% [epoch badChan]=bciGetEpoch(bci,dat)
% creates a single epoch of data segment according to parameter setting in
% bciGetParam (which is recommended to be located in the subject's folder)

function [epoch badChan]=bciGetEpoch(bci,dat)

% 110502 CC included bci.init.nBaseline>1 otherwise code not executed,
%        implemented Projected common average reference

% cast double precision, otherwise absurd results with filtfilt and maybe
% others
dat = double(dat);

% do baseline correction
if bci.init.nBaseline>=1
    baseline = repmat(mean(dat(:,1:bci.init.nBaseline),2),1,bci.init.nEpochLength);
    dat = dat - baseline;
end


% detrend time series
if bci.param.detrend,
    dat=detrend(dat')';
end

% find artifacts
badChan = sum(abs(dat(:,bci.init.intervalOfInterest))>bci.param.artifactThresh,2) > 0;

% remove common average reference
if ~isempty(bci.refChanCAR)
    ecog.data = dat;
    % ecog.selectedChannels = bci.preprocChan; % preprocChan already
    % applied in call of this function
    ecog.selectedChannels = 1:size(dat,1); % exclude bad channels of this single epoch???
    %ecog.selectedChannels = setdiff(1:size(dat,1),find(badChan));
    ecog=ecogRemoveCommonAverageReference(ecog,'projected');
    
    dat = ecog.data;
%     refChan=setdiff(bci.refChanCAR,find(badChan));
%     if isempty(refChan),
%         refChan = bci.refChanCAR;
%     end
%     dat = dat - ones(size(dat,1),1)*mean(dat(refChan,:),1);
end

% spatial filter
if bci.param.spatialFilter>0,
    dat = (dat'*bci.init.spatialFilter)';
end

% bandpass filter, important: cast to double 
if length(bci.param.bandFreqs)>0,
    dat = filtfilt(bci.init.filtCoeffB,bci.init.filtCoeffA,dat')';
end

if bci.param.specMethod==1,
    % Fourier transform
    epoch = bciGetEpochFFT(bci,dat);
elseif  bci.param.specMethod==2,
    % multitaper
    tmp = mtspectrumc(dat(:,bci.init.intervalOfInterest)',bci.init.taperParams)';
    epoch=tmp(:,bci.init.fftFreqIdx);
    % temporal multitaper?
elseif  bci.param.specMethod==3,
    % AR model
   tmp= lpc(dat(:,bci.init.intervalOfInterest)',bci.param.ARorder);
   epoch = tmp(:,2:end);
elseif  bci.param.specMethod==0,
    % timeseries data
    if ~isempty(bci.param.resampFreq),
        % resample
        dat=resample(dat', floor(bci.param.resampFreq), floor(bci.srate))';
        intOfI=bci.init.resampIntOfInterest;
    else
        intOfI=bci.init.intervalOfInterest;
    end
    epoch=dat(:,intOfI);
elseif  bci.param.specMethod==4,
    % exact FT
    epoch = bciGetEpochExactFT(bci,dat);
elseif  bci.param.specMethod==5,
    % exact Taper
    epoch = bciGetEpochExactTaper(bci,dat);
elseif  bci.param.specMethod==6,
    % fft based SNR (development status)
    epoch = bciGetEpochSNRfft(bci,dat);
elseif  bci.param.specMethod==7,
    % exact FT based SNR (development status)
    epoch = bciGetEpochSNRexact(bci,dat);
elseif bci.param.specMethod==8,
    % combine exact fft and exact SNR (development status)
    epoch = cat(2,bciGetEpochExactFT(bci,dat),bciGetEpochSNRexact(bci,dat));
elseif bci.param.specMethod==9,
    % combine exact taper and exact SNR (development status)
    epoch = cat(2,bciGetEpochExactTaper(bci,dat),bciGetEpochSNRexact(bci,dat));
elseif bci.param.specMethod==10,
    % combine fft and exact taper (development status)
    epoch = cat(2,bciGetEpochFFT(bci,dat,1:4),bciGetEpochExactTaper(bci,dat,5:8));
elseif bci.param.specMethod==11,
    % combine exact taper and exact taper (development status)
    epoch = cat(2,bciGetEpochExactTaper(bci,dat,1:4),bciGetEpochExactTaper(bci,dat,5:8));    
else
    error('bci.param.specMethod not defined');
end

%------------------------------
%% feature extraction functions
%------------------------------
function epoch = bciGetEpochFFT(bci,dat, destFreq)
% get fft data
Y = fft(dat(:,bci.init.intervalOfInterest).*(ones(size(dat,1),1)*bci.init.fftWin),bci.init.NFFT,2);
tmp=abs(Y).^2;
if bci.param.fftSmooth,
    average filter
    tmpConv=conv2(tmp,[1 1 1]/3);
    tmp=tmpConv(:,2:end-1);
end
if nargin ==2,
    epoch=tmp(:,bci.init.fftFreqIdx);
else
    epoch=tmp(:,bci.init.fftFreqIdx(destFreq)); % internal frequency selection (for combined feature extraction)
end

%--------------------------------------------------------------------------

function epoch = bciGetEpochExactFT(bci,dat)
dataProj = repmat(dat(:,bci.init.intervalOfInterest),[1,1,length(bci.param.freq)]);
epoch=squeeze(abs(mean(dataProj.*bci.init.eFT.eFunc(1:size(dat,1),:,:,1),2)));

%--------------------------------------------------------------------------
    
function epoch = bciGetEpochExactTaper(bci,dat, destFreq)
if nargin ==2,
    destFreq=1:length(bci.param.freq);
end
if size(bci.param.tapers,1)>1,
    tap = unique(bci.param.tapers(destFreq,2));
else
    tap = size(bci.init.eFT.eFunc,4);
end
dataProj = dat(:,bci.init.intervalOfInterest,ones(tap,1));
dataProj = repmat(permute(dataProj,[1 2 4 3]),[1 1 length(destFreq),1]);
if 0, % absolute value
    epoch=abs(mean(dataProj.*bci.init.eFT.eFunc(1:size(dat,1),:,destFreq,1:tap),2));
    epoch = squeeze(mean(epoch,4));
else % complex conjugate
    epoch=mean(dataProj.*bci.init.eFT.eFunc(1:size(dat,1),:,destFreq,1:tap),2);
    epoch = squeeze(mean(conj(epoch).*epoch,4));
end

%--------------------------------------------------------------------------

function epoch = bciGetEpochSNRfft(bci,dat)
% SNR of fourier transform
% 1st approach: fast fourier transform
Y = fft(dat(:,bci.init.intervalOfInterest).*(ones(size(dat,1),1)*bci.init.fftWin),bci.init.NFFT,2);
tmp=abs(Y).^2;
noiseWidth = 2.0; % bandwidth of noise in Hz
nSteps = max(1,round(noiseWidth / (bci.srate/bci.init.NFFT)));
noise = zeros(size(tmp,1),length(bci.init.fftFreqIdx));
for k=1:nSteps,
    noise = noise + tmp(:,bci.init.fftFreqIdx+k) + tmp(:,bci.init.fftFreqIdx-k);
end
epoch=tmp(:,bci.init.fftFreqIdx)./(noise/2/nSteps);

%--------------------------------------------------------------------------

function epoch = bciGetEpochSNRexact(bci,dat)
% 2nd approach: exact fourier transform
dataProj = repmat(dat(:,bci.init.intervalOfInterest),[1,1,length(bci.param.freq)]);
epoch=squeeze(abs(mean(dataProj.*bci.init.eFT.eFunc(1:size(dat,1),:,:,1),2)));
noise = zeros(size(epoch));
noiseSteps = [0.2:0.2:2.0];
for fi = 1:size(epoch,2),
    idx = 1:bci.init.eFT.periodMaxIdx(fi);
    t=linspace(0,bci.param.winsize,bci.init.nWinLength);
    for k=1:length(noiseSteps),
        noise(:,fi) = noise(:,fi) +  ...
            abs(mean(dataProj(:,idx,fi).*(ones(size(dat,1),1)*exp(-1i*2*pi*(bci.param.freq(fi)+noiseSteps(k))*t(idx))),2)) + ...
            abs(mean(dataProj(:,idx,fi).*(ones(size(dat,1),1)*exp(-1i*2*pi*(bci.param.freq(fi)-noiseSteps(k))*t(idx))),2));
    end
end
epoch = epoch./(noise/2/length(noiseSteps));

        